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ABSTRACT 


We  present  a  straightforward  methodology  for  converting  the  deterministic  multi-wave  pa- 
rameterizations  of  nonorographic  gravity-wave  drag,  currently  used  in  general  circulation 
models  (GCMs),  to  stochastic  analogues  that  use  far  fewer  waves  (in  our  example,  a  single 
wave)  within  each  grid  box.  Deterministic  discretizations  of  source-level  momentum  flux 
spectra  using  a  fixed  spectrum  of  many  waves  with  predefined  phase  speeds  are  replaced 
by  sampling  these  source  spectra  stochastically  using  waves  with  randomly-assigned  phase 
speeds.  Using  simple  conversion  formulas,  we  show  that  time-mean  wave-induced  drag, 
diffusion  and  heating-rate  profiles  identical  to  those  from  the  deterministic  scheme  are  pro¬ 
duced  by  the  stochastic  analogue.  Furthermore,  the  need  for  bulk  intermittency  factors  of 
small  value  is  largely  obviated  through  the  explicit  incorporation  of  stochastic  intermittency 
into  the  scheme.  When  implemented  in  a  GCM,  the  single-wave  stochastic  analogue  of 
an  existing  deterministic  scheme  reproduces  almost  identical  time-mean  middle-atmosphere 
climate  and  drag  as  its  deterministic  antecedent,  but  with  an  order  of  magnitude  reduc¬ 
tion  in  computational  expense.  The  stochastic  parameterization  is  accompanied  by  natural 
stochastic  variability  about  the  time-mean  profile  that  forces  the  smallest  space-time  scales 
of  the  GCM.  Studies  of  mean  GCM  kinetic  energy  spectra  show  that  this  additional  stochas¬ 
tic  forcing  does  not  lead  to  unrealistic  increases  in  dynamical  variability  at  these  smallest 
GCM  scales,  although  a  systematic  increase  in  divergent  kinetic  energy  variability  is  evident. 
Our  results  show  that  the  expensive  deterministic  schemes  currently  used  in  GCMs  are  easily 
modified  and  replaced  by  cheap  stochastic  analogues  without  any  obvious  deliterious  impacts 
on  GCM  climate  or  variability,  while  offering  potential  advantages  of  computational  savings, 
reduction  of  systematic  climate  biases  and  greater  and  more  realistic  ensemble  spread. 
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1.  Introduction 


Finite  computational  resources  force  global  weather  and  climate  prediction  models  to 
run  at  spatial  resolutions  that  do  not  resolve  the  full  spectrum  of  gravity  waves  that  can 
exist  in  the  atmosphere.  Since  the  dissipation  of  gravity-wave  momentum  and  energy  induce 
significant  body  forces,  heating  and  constituent  mixing  at  synoptic  scales,  general  circulation 
models  (GCMs)  must  parameterize  these  missing  gravity- wave- induced  effects  on  the  resolved 
flow  (Kim  et  al.  2003).  Parameterizations  of  drag  due  to  unresolved  orographic  gravity  waves 
were  first  implemented  in  weather  and  climate  models  over  two  decades  ago,  where  they 
had  immediate  positive  influences  in  the  winter  extratropical  troposphere  and  stratosphere 
(Palmer  et  ah  1986;  McFarlane  1987).  They  are  now  essential  components  of  any  credible 
global  weather  or  climate  prediction  system. 

Parameterizations  of  gravity  waves  from  nonorographic  sources  were  longer  in  coming, 
despite  emerging  understanding  of  their  primary  role  in  controlling  the  large-scale  circula¬ 
tion  of  the  middle  atmosphere,  particularly  in  the  tropics  and  summer  extratropics  (Lindzen 
and  Holton  1968;  Dunkerton  1982b;  Holton  1983;  Garcia  and  Solomon  1985).  Develop¬ 
ment  was  stymied  at  first  by  insufficient  observational  and  theoretical  knowledge  of  relevant 
nonorographic  wave  sources.  High-resolution  observations  of  gravity  wave-induced  velocity 
and  temperature  perturbations  later  revealed  a  broad  spectrum  of  waves  throughout  the 
troposphere  and  middle  atmosphere  with  surprisingly  reproducible  spectral  shapes  (Smith 
et  ah  1987),  which  motivated  an  initial  generation  of  nonorographic  gravity  wave  drag  pa¬ 
rameterizations  based  on  a  quasi- invariant  global  background  spectrum  of  many  waves  from 
indistinct  tropospheric  sources  (e.g.,  Fritts  and  VanZandt  1993;  Warner  and  McIntyre  1996). 
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A  variety  of  these  spectral  nonorographic  gravity- wave  drag  schemes  now  exist  (e.g.,  Kiehl 
et  ah  1996;  Alexander  and  Dunkerton  1999;  Medvedev  and  Klaassen  2000;  Garcia  et  al.  2007) 
and  they  constitute  the  standard  means  of  parameterizing  nonorographic  gravity  wave  drag 
in  global  models  at  present  (see,  e.g.,  Table  1  of  Eyring  et  al.  2006).  A  next  generation 
of  schemes  is  slowly  emerging,  based  on  physical  models  of  gravity-wave  generation  from 
specific  nonorographic  sources  such  as  deep  convection  and  frontogenesis:  they  too  launch 
a  broad  spectrum  of  gravity  waves  (e.g.,  Charron  and  Manzini  2002;  Song  and  Chun  2005) 
and  can  thus  often  be  implemented  by  simply  replacing  the  uniform  source-level  momentum 
flux  function  of  a  pre-existing  spectral  nonorographic  scheme  (e.g.,  Beres  et  al.  2005;  Richter 
et  al.  2010). 

These  parameterizations  of  nonorographic  gravity-wave  drag  typically  specify  source- 
level  wave  momentum  flux  as  a  function  of  ground-based  horizontal  phase  speed  c,  denoted 
Arc(c),  which  is  then  discretized  among  ngw  individual  gravity  waves  of  phase  speed  Cj  and 
momentum  flux  Tj ,  where  j  —  1 .  .  .ngw.  After  assigning  the  remaining  parameters  of  each 
tagged  wave  j  (e.g.,  horizontal  wavenumber  vectors  Kj),  the  propagation  and  dissipation 
modules  then  determine  how  each  wave’s  momentum  flux  is  deposited  into  higher  model 
levels.  The  resulting  tendencies  due  to  all  the  waves  are  summed  and  then  applied  to  modify 
model  winds  and  temperatures. 

These  rsrc(c)  functions  are  typically  broad,  so  that  a  large  number  of  waves  ngw  is  often 
required  for  a  sufficiently  accurate  discretization.  Thus,  unlike  orographic  gravity  wave 
parameterizations  that  typically  launch  only  1  or  2  waves  in  each  model  grid  box  (Scinocca 
and  McFarlane  2000;  Webster  et  al.  2003),  nonorographic  schemes  can  launch  anywhere  from 
ngw  ~  1 0  1000  parameterized  waves  (Alexander  and  Dunkerton  1999;  Scinocca  2003;  Garcia 
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et  al.  2007;  Zhu  et  al.  2010;  Orr  et  al.  2010).  Consequently,  nonorographic  gravity  wave 
drag  parameterizations  can  be  computationally  expensive,  which  has  spurred  recent  efforts 
to  speed  up  specific  schemes  to  make  practical  their  integration  into  general  circulation 
models  (GCMs)  used  for  production  runs  (e.g.,  Warner  and  McIntyre  2001;  Scinocca  2003). 
This  generally  involves  a  set  of  simplifications  or  optimizations  specific  to  that  particular 
parameterization  that  do  not  change  the  underlying  algorithm  or  output  in  any  major  way. 

Here  we  investigate  a  different  approach  to  this  issue  that  is  potentially  applicable  to 
any  existing  multiwave  parameterization  of  nonorographic  gravity-wave  drag.  The  central 
idea  is  to  replace  the  deterministic  discretizations  of  rsrc(c)  into  ngw  individual  waves  with 
a  stochastic  discretization  that  can  involve  just  a  single  parameterized  wave  in  each  model 
grid  box.  The  approach  is  developed  mathematically  in  application  to  a  specific  multiwave 
parameterization  of  nonorographic  gravity- wave  drag  described  in  section  2.  The  stochastic 
analogue  is  described  in  section  3  and  compared  to  its  deterministic  parent  in  single-column 
tests  in  section  4.  The  two  are  implemented  in  a  GCM  in  section  5  and  the  GCM  climate 
and  variability  that  result  from  each  in  long-term  integrations  are  compared  and  contrasted. 
The  results  are  discussed  in  section  6  and  the  major  findings  and  implications  summarized 
in  section  7. 


2.  Deterministic  Parameterization 

While  the  ideas  to  follow  are  general,  we  illustrate  and  implement  them  here  for  one 
specific  scheme:  the  multiwave  parameterization  of  nonorographic  gravity  wave  drag  imple¬ 
mented  in  version  3.0  of  the  Whole  Atmosphere  Community  Climate  Model  (WACOM),  as 
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summarized  in  Appendix  A  of  Garcia  et  al.  (2007).  Full  details  of  the  scheme’s  formulation 
and  numerics  are  provided  by  Kichl  et  al.  (1996),  Collins  et  al.  (2004)  and  Garcia  et  al. 
(2007),  while  Appendix  A  outlines  the  recent  generalization  of  this  parameterization  code 
into  a  “team  scheme”  for  use  in  GCMs  at  various  US  institutions.  Here  we  discuss  only 
those  aspects  of  the  scheme  salient  to  the  present  work. 

The  scheme’s  prescribed  source-level  momentum  flux  function  takes  the  Gaussian  form 


Tsrc(c)  =  n  exp 


(C  Coff ) 


Tb  =  TbF{(j>,t), 


(1) 

(2) 


with  a  phase-speed  width  cw  =  30  m  s^1.  rb  is  the  “background”  momentum  flux  and  is 
scaled  from  the  constant  baseline  value  r6*  by  a  factor  F((j),t),  which  varies  with  latitude 
cf)  and  time  (season)  t,  as  given  by  eqs.  (A24)-(A26)  of  Garcia  et  al.  (2007)  and  plotted  in 
Figure  la. 

This  source  flux  rsrc(c)  is  discretized  at  the  launch  pressure  level  psrc  =  500  hPa  by 
assigning  an  equispaced  distribution  of  ground-based  horizontal  phase  speeds 


Cj  ^off  T  j*Ac,  (3) 

j*  =  ~nc,  —nc  +  1, . . .  +  nc  -  1,  +nc  =  j  -  nc  +  1,  (j*  e  Z,  nc  G  N) , 

where  Ac  is  the  phase  speed  resolution,  yielding  a  total  of  ngw  =  2 nc  +  1  individual  gravity 
waves  of  momentum  flux  Tsrc(cj).  This  phase  speed  distribution  is  symmetric  about  caff  and 
thus  samples  (1)  symmetrically  about  it’s  peak. 

In  the  Garcia  et  al.  (2007)  formulation,  c0ff  =  Usrc  =  Usrc| ,  where  Usrc  is  the  horizontal 
velocity  vector  at  psrc.  Horizontal  wavenumber  vectors  Ky  align  parallel  to  Usrc,  yielding 
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a  symmetric  distribution  of  intrinsic  horizontal  phase  speeds  with  flux  peaking  near  zero 
intrinsic  phase  speed.  Garcia  et  al.  (2007)  assign  ngw  =  65  waves  with  Ac  =  2.5  m  s_1, 
which  yields  intrinsic  phase  speeds  spanning  the  range  ±80  m  s-1.  The  resulting  discretized 
sampling  of  the  normalized  flux  function  Tsrc(c)/Tb  is  plotted  in  black  in  Figure  2a.  For 
comparison,  the  ngw  =  9,  Ac  =  10  nr  s^1  discretization  used  by  Kiehl  et  al.  (1996)  is  shown 
in  gray,  which,  since  they  set  cQ//  =  0,  spans  a  ±40  m  s_1  range  of  ground-based  phase 
speeds. 

The  subsequent  deposition  of  wave  momentum  flux  at  higher  altitudes  is  parameterized 
for  each  wave  using  a  Lindzen  (1981)  parameterization  of  hydrostatic  irrotational  vertical 
propagation  subject  to  critical-level  removal  and  linear  saturation  thresholds  (see  Kiehl  et  al. 
1996;  Garcia  et  al.  2007,  for  details).  The  ensuing  mean-flow  acceleration  (or  gravity 
wave  drag  per  unit  mass)  at  model  layer  k  due  to  wave  j  of  phase  speed  Cj  is 

<9r,-  k 

aj,k  = 

op 

where  g  is  gravitational  acceleration,  p  is  pressure,  and  e  is  a  constant  in  the  range  0  <  e  <  1 
that  represents  the  “intermittency”  or  “efficiency”  of  wave  breaking.  Intermittency  factors 
of  this  sort  appear  in  many  multiwave  nonorographic  gravity-wave  drag  schemes  but  their 
implementation  and  effects  can  vary  from  scheme  to  scheme.  The  implementation  here 
follows  that  in  Alexander  and  Dunkerton  (1999),  and  so,  as  in  their  scheme,  if  we  retain  the 
same  ±80  m  s-1  range  of  phase  speeds,  then  e  scales  linearly  with  changes  in  Ac  and  inversely 
with  the  corresponding  changes  in  ngw  to  retain  the  same  total  mean-flow  acceleration. 

That  total  mean-flow  acceleration 

‘klgw 

ak  =  aj,ki  (5) 

1=1 
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directed  parallel  to  Usrc,  is  projected  into  zonal  and  meridional  components  that  are  applied 
as  tendencies  to  the  model’s  horizontal  velocity  field  U*,.  The  WACOM  3.0  scheme  imple¬ 
ments  as  options  a  number  of  additional  limits  on  single-wave  and  total  tendencies.  These 
limits  are  applied  in  such  a  way  that  flux  is  redistributed  rather  than  removed/added,  so 
as  to  conserve  column- integrated  momentum.  The  version  used  here  deposits  all  remaining 
wave  flux  in  the  top  two  model  layers  to  ensure  robust  circulation  and  climate  responses 
(Shaw  et  al.  2009). 

With  accelerations  specified,  other  quantities  follow  based  on  the  Lindzen  (1981)  satu¬ 
ration  model.  The  effective  vertical  diffusion  coefficient  due  to  the  turbulence  generated  by 
wave  breaking  is 

Pgw 

{Dgwd)k  =  yy 

j= i 

where  Upro:>  is  the  component  of  the  wind  vector  U*,  projected  along  Kj,  Nk  is  buoyancy 
frequency  in  layer  k,  and  Pr  is  the  effective  Prandtl  number,  here  set  equal  to  4  following 
Garcia  et  al.  (2007). 

The  wave- induced  heating  rate  employed  here  is 


Pr 


-n 


-ur  ij 


2  Nl 


O'jik 


(6) 


m 


ngw  s 

E  c-  k  ~  ur) 

3= 1  1  P 


d 


dt  ^  1  C„  V~J  y  aj,k  (1  +  Pr)  dp 

and  is  based  on  the  work  of  Medvedev  and  Klaassen  (2003),  where  Tk  and  pk  are  temperature 
and  density,  respectively,  in  layer  k  and  Cp  is  mass  specific  heat  at  constant  pressure.  The  first 
term  is  a  uniformly  positive  irreversible  heating  term  due  to  deposition  of  total  wave  energy, 
both  the  frictional  dissipation  of  wave  kinetic  energy  and  the  thermal  dissipation  of  wave 
potential  energy.  The  second  is  a  differential  heating/cooling  term  associated  with  vertical 
variations  in  the  wave  heat  flux,  which  Akmaev  (2007)  shows  can  only  result  from  thermal 


[Pk  fe  -  ur)  «**]  , 


(7) 


7 


dissipation  of  wave  potential  energy,  leading  to  the  1/(1  +  Pr)  factor  in  (7).  This  heating 
rate  expression  is  the  only  part  of  the  parameterization  used  here  that  differs  substantially 
from  that  described  by  Garcia  et  ah  (2007). 


3.  Stochastic  Analogue 

Absent  a  specific  nonorographic  source  model,  the  “background”  flux  function  (1)  is 
a  practical  choice  that  simplifies  fitting  of  source  parameters  to  observed  climatological 
distributions  of  gravity- wave  phase  speeds  and  momentum  fluxes  (see,  e.g.,  Alexander  and 
Vincent  2000;  Gong  et  ah  2008).  To  motivate  what  follows,  however,  a  physical  interpretation 
of  it  is  useful  here.  For  present  purposes,  one  can  view  it  as  the  state  that  emerges  once  a 
large  number  of  random  wave  modes  attain  some  form  of  statistical  mechanical  equilibrium 
(an  analogy  pursued  explicitly  in  some  spectral  gravity- wave  models:  e.g.,  Allen  and  Joseph 
1989;  Souprayen  et  ah  2001).  The  equilibrium  spectrum  (1)  would  then  emerge  only  over  a 
volume  and  time  both  large  and  long  enough,  respectively,  for  the  full  ensemble  of  gravity- 
wave  wavelengths  and  periods  to  attain  equilibrium. 

Given  gravity-wave  horizontal  wavelengths  of  up  to  1000  km  and  periods  and  group- 
propagation  times  of  up  to  a  day,  typical  GCM  grid-box  dimensions  of  10-1000  km  and 
time  steps  of  1-60  min  would  not  appear  to  be  either  large  or  long  enough,  respectively, 
for  this  wave  ensemble  to  equilibrate  spectrally.  Spectral  equilibrium  would  appear  instead 
only  when  averaged  over  wider  horizontal  areas  encompassing  many  GCM  grid  boxes,  and 
when  averaged  over  a  number  of  GCM  time  steps.  At  any  given  time  step  in  one  grid  box, 
subgridscale  wave  flux  would  instead  resemble  the  random  nature  of  the  component  wave 


modes  themselves.  Clearly,  depending  on  model  resolution,  there  is  a  continuum  of  possible 


states  within  a  grid  box,  with  the  background  multiwave  “equilibrium”  spectrum  and  a 
purely  random  stochastic  state  representing  the  two  end  limits. 

A  very  simple  approach  to  parameterizing  such  states  is  depicted  in  Figure  2b.  Instead 
of  discretizing  (1)  deterministically  with  ngw  equispaced  wave  phase  speeds  (Figure  2a),  we 
now  sample  it  randomly  by  choosing  nsgw  “stochastic”  waves  with  phase  speeds 


cQff  ~t~  cr  (2 Rj  1) , 


(8) 


where  j  =  1 . .  .nsgw,  cr  =  80  m  s'1  is  the  phase  speed  range  and  Rj  is  the  output  from 
a  random  number  generator  with  a  uniform  mean  distribution,  such  that  0  <  Rj  <  1.  In 
this  implementation,  the  random  Rj  values  are  repopulated  at  every  grid  point  and  at  every 
model  time  step  so  that  there  are  no  spatiotemporal  correlations  in  wave  properties  between 
adjacent  grid  boxes  or  model  time  steps,  in  contrast  to  the  original  deterministic  scheme  in 
which  waves  at  adjacent  grid  boxes  and  times  are  highly  correlated. 

For  a  given  bulk  intermittency  e  used  in  the  original  deterministic  scheme,  the  same  time 
mean  momentum  fluxes  and  total  accelerations  are  attained  in  the  stochastic  analogue  by 
using  a  scaled  intermittency  in  (4)  of 


(9) 


If  e  parameterizes  the  bulk  effects  of  stochastic  intermittency  alone  in  the  deterministic 
scheme,  one  could  argue  such  factors  should  be  removed  entirely  from  the  stochastic  analogue 
in  which  this  intermittency  is  now  explicit.  Doing  so  is  straightforward  by  adding  a  second 


uniformly  distributed  random  variable  Sj  (0  <  Sj  <  1)  and  choosing  a  limit  S ,  such  that  if 
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Sj  <  S  the  stochastic  acceleration  aJ(ot  is  applied  in  the  model,  but  if  Sj  >  S  we  set  a j)ot  =  0. 
Then  es  disappears  from  (4),  and  is  replaced  by  the  new  stochastic  analogue  of  (9): 

S  =  (10) 

X^sgw  ) 

Note  that  eqs.  (9)  and  (10)  can  yield  values  in  excess  of  unity.  While  es  >  1  is  technically 
unphysical  the  stochastic  parameterization  algorithm  still  works  using  such  es  settings.  By 
contrast  S  >  1  cannot  be  accomodated,  and  so  nsgw  needs  to  be  increased  until  S  <  1  is 
achieved.  Another  way  of  viewing  (10)  is  as  a  generalization  of  the  stochastic  scheme  to  a 
noninteger  number  of  waves 


n 


sgw 


^sgw^  • 


(11) 


This  permits,  for  example,  implementations  with  less  than  one  stochastic  wave  per  gridbox 
(0  <  hsgw  <  1),  by  choosing  nsgw  =  1  and  nonvanishing  S  <  1.  For  simplicity,  in  this  paper 
we  only  show  results  using  the  es  formulation,  given  the  more  straightforward  connection  to 
its  deterministic  antecedent. 


4.  Offline  Single  Column  Tests 

A  convenient  feature  of  the  stochastic  implementation  in  section  3  is  the  close  connection 
that  is  maintained  to  the  original  deterministic  scheme.  The  latter  has  been  carefully  refined 
and  tuned  for  use  in  global  models  over  many  years.  The  simple  relations  in  section  3  allow 
the  core  physics  and  tuned  parameter  settings  of  the  determinstic  scheme  to  translate  to  the 
stochastic  analogue,  which  should  in  turn  yield  the  same  mean  drag  profiles.  We  demonstrate 
this  here  using  offline  single-column  tests. 
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The  top  panels  of  Figure  3  show  vertical  profiles  of  instantaneous  zonal  and  meridional 
winds  at  a  grid  point  near  the  Alps  after  +12  hours  of  a  T79L68  global  model  forecast 
initialized  on  1  June  2007  at  0000  UTC.  The  model  in  question  is  the  Advanced-Level 
Physics  High- Altitude  prototype  of  the  Navy  Operational  Global  Atmospheric  Prediction 
System  (NOGAPS- ALPHA),  and  the  fields  are  part  of  the  high- altitude  forecast-assimilation 
runs  described  by  Eckermann  et  al.  (2009). 

The  forecast  model  in  those  runs  used  the  deterministic  WACCM  3.0  scheme  to  param¬ 
eterize  nonorographic  gravity  wave  drag.  Thus  at  this  time  and  location  the  model  passed 
these  exact  wind  (and  other  meteorological)  profiles  to  that  parameterization,  which  in  turn 
returned  zonal  and  meridional  mean-flow  accelerations,  vertical  diffusivities  and  dynamical 
heating  rates  shown  with  thick  solid  gray  curves  in  the  remaining  panels  of  Figure  3.  Here 
we  have  used  the  same  tuned  parameter  settings  as  in  Eckermann  et  al.  (2009),  specifically 
ngw  =  65,  Ac  =  2.5  nr  s-1,  psrc  =  500  hPa,  rb*  =  1.75  nrPa  and  e  =  0.0175. 

Other  curves  in  these  lower  four  panels  show  results  from  our  stochastic  analogue  of  this 
scheme  that  uses  only  a  single  wave  ( nsgw  =  1).  Here  output  from  the  stochastic  parame¬ 
terization  was  averaged  over  a  number  of  separate  calls  ranging  from  1  to  10000,  using  the 
same  input  wind  profiles  in  every  case,  but  with  the  random  number  Rj  in  (8)  independently 
reinitialized  during  each  call.  The  mean  accelerations,  diffusivities  and  heating  rates  after 
1-10  calls  in  Figures  3c-3f  are  substantially  different  from  the  deterministic  reference  (gray 
curve)  due  to  the  random  nature  of  the  wave  field.  However,  after  100  calls  the  mean  profiles 
are  quite  similar  to  the  reference  curve,  and  after  1000-10000  calls  the  mean  profiles  over¬ 
lay  the  deterministic  reference  curves.  However,  as  shown  in  Figure  4,  while  the  long-term 
means  are  the  same,  the  stochastic  version  produces  large  standard  deviations  about  that 
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mean,  whereas,  for  the  same  input  meteorological  profiles,  the  original  deterministic  scheme 
has  zero  standard  deviation. 

Since  nsgw  =  1  and  ngw  =  65,  from  (9)  we  used  a  modified  stochastic  intermittency 
es  =  65e  ~  1.138.  Thus,  by  incorporating  intermittency  into  the  parameterization  explicitly 
(Figure  4),  the  need  for  a  parameterized  bulk  intermittency  factor  is  now  largely  obviated. 
Of  course  the  es  ~  1  result  is  specific  to  the  tuned  settings  for  this  particular  model  config¬ 
uration  and  thus  probably  fortuitous,  with  tuned  values  for  other  models  likely  leading  to 
es  7^  1.  Nonetheless,  the  trend  away  from  very  small  e  values,  implying  highly  intermittent 
or  inefficient  wave  breaking,  to  values  nearer  unity  through  an  explicit  stochastic  repre¬ 
sentation  of  intermittency  in  the  parameterization,  is  clearly  both  a  robust  and  physically 
self-consistent  result. 

One  can  reduce  the  large  standard  deviations  in  Figure  4  by  increasing  nsgw.  Figure  5 
shows  the  corresponding  mean  zonal  accelerations  and  standard  deviations  as  nsgw  is  pro- 
gessively  increased,  with  es  rescaled  in  each  case  as  in  (9).  As  nsgw  increases,  the  standard 
deviation  reduces  towards  the  vanishing  deterministic  limit.  Thus  nsgw  >  1  yields  hybrid 
states  that  are  a  blend  of  the  stochastic  ( nsgw  =  1)  and  deterministic  ( nsgw  — >  oo)  limits. 


5.  Global  Model  Tests 

Next  we  compare  how  the  equivalent  stochastic  and  determinstic  versions  of  this  nonoro- 
graphic  gravity-wave  drag  parameterization  perform  in  a  GCM.  We  use  the  forecast  model 
component  of  NOGAPS-ALPHA  with  the  same  T79L68  formulation  and  physics  settings 
described  by  Eckermann  et  ah  (2009),  except  that  here  we: 


12 


•  use  a  Webster  et  al.  (2003)  parameterization  of  orographic  gravity-wave  and  flow¬ 
blocking  drag  instead  of  the  Palmer  et  al.  (1986)  scheme. 

•  apply  the  nonorographic  gravity  wave-induced  heating  rate  (7)  to  GCM  temperature 
fields  using  Pr  =  4. 

The  GCM  is  intialized  on  1  June  2007  using  the  NOGAPS- ALPHA  analysis  fields  described 
by  Eckermann  et  al.  (2009),  then  is  integrated  forward  in  time  without  assimilation  update 
cycles  to  1  January  2010.  This  “nature  run”  is  constrained  by  12  hourly  analyzed  sea-  and 
land-surface  temperatures,  snow  depths  and  ice  concentrations  at  the  lower  boundary. 

A  series  of  these  nature  runs  was  performed  initially  to  tune  the  nonorographic  gravity- 
wave  drag  parameterization  to  yield  a  realistic  zonal-mean  middle  atmosphere  climate,  which 
led  to  several  changes  from  the  default  WACCM  3.0  settings  described  in  earlier  sections. 
First,  the  background  flux  t&  in  (2)  was  modified,  as  shown  in  Figure  lb,  to  center  flux 
peaks  closer  to  the  solstices  and  to  increase  wave  momentum  fluxes  at  the  equator.  Second, 
we  chose  to  launch  nonorographic  waves  zonally  rather  than  along  the  source-level  wind 
direction.  Third,  we  reduced  the  critical  inverse  Froude  number  for  nonorographic  gravity- 
wave  breaking,  Fr~x,  from  1  to  0.1,  to  force  parameterized  waves  to  break  at  lower  altitudes, 
an  approach  often  used  to  tune  both  orographic  and  nonorographic  gravity-wave  drag  in 
GCMs  (e.g.,  Norton  and  Thuburn  1999;  Webster  et  al.  2003;  Scinocca  et  al.  2008)  and 
recently  defended  on  theoretical  grounds  by  Scinocca  and  Sutherland  (2010). 

Figure  6a  shows  3-year  average  zonal-mean  zonal  winds  for  July  from  a  control  run  with¬ 
out  parameterized  nonorographic  gravity-wave  drag,  revealing  unrealistically  strong  strato¬ 
spheric  jets  that  extend  through  the  mesosphere.  Figure  6b  shows  corresponding  mean  winds 


13 


from  the  nature  run  in  which  the  deterministic  nonorographic  gravity-wave  drag  was  acti¬ 
vated  using  the  tuned  settings  noted  above,  with  rfe*  =  10  mPa,  e  =  0.0375  and  ngw  =  65. 
These  simulations  with  tuned  nonorographic  gravity-wave  drag  show  more  realistic  strato¬ 
spheric  jets  in  both  hemispheres,  including  better  tilting  of  the  winter  (southern)  jet  and 
realistic  reversal  of  the  summer  (northern)  jet  in  the  upper  mesosphere. 

Figure  6c  shows  results  from  the  nature  run  using  the  stochastic  analogue  of  the  tuned 
deterministic  nonorographic  gravity-wave  drag,  using  nsgw  =  1  and  thus  es  =  2.275  according 
to  (9).  Visual  comparison  of  Figures  6b  and  6c  reveals  almost  identical  zonal  wind  structure, 
despite  the  imposition  of  explicitly  stochastic,  highly  intermittent  and  “noisy”  gravity-wave 
drag  and  heating  rates  in  the  latter  GCM  simulation.  To  verify  this  visual  impression, 
Figure  6d  plots  the  mean  zonal  wind  difference  fields  between  the  stochastic  and  deterministic 
nature  runs.  Differences  everywhere  are  small,  particularly  in  the  summer  hemisphere. 

Figure  7  plots  the  3-year  zonal-mean  zonal  mean-flow  accelerations  (top  row)  and  heating 
rates  (bottom  row)  for  July  due  to  parameterized  nonorographic  gravity- wave  drag  from  the 
deterministic  and  stochastic  nature  runs.  The  time-mean  accelerations  and  heating  rates  in 
the  GCM  again  look  largely  identical.  Small  values  of  the  difference  fields,  plotted  on  the 
right  of  Figure  7,  again  confirm  that  impression. 

Since  our  offline  single  column  simulations  demonstrated  that  the  stochastic  approach 
gave  identical  time-mean  accelerations  and  heating  rates  to  its  determinstic  antecedent,  these 
GCM  results  might  not  seem  all  that  surprising.  Yet  in  a  fully  interactive  nonlinear  GCM, 
it  is  not  a  given  that  highly  intermittent  stochastic  drag  will  produce  the  same  long-term 
GCM  climate  as  it’s  smooth  deterministic  progenitor.  Indeed,  the  corresponding  two-year 
zonal-mean  zonal  winds  in  January,  plotted  in  Figure  8,  show  that  the  winter  (northern) 
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stratospheric  zonal  winds  in  this  case  are  very  different  between  the  stochastic  and  deter¬ 
ministic  nature  runs.  These  differences  arose  due  to  spontaneous  stratospheric  warmings  in 
both  Januaries  of  the  stochastic  run  which  did  not  occur  in  the  deterministic  simulation:  the 
latter  generated  a  warming  only  in  December  2007,  which  did  not  occur  in  the  stochastic 
run.  Of  course,  stratospheric  warmings  in  northern  winter  are  a  well-known  source  of  natural 
interannual  variability  in  GCMs,  and  nature  runs  extending  for  25-50  years  or  more  would 
be  needed  to  deduce  any  real  systematic  differences  in  zonal-mean  northern  winter  strato¬ 
spheric  climate  or  stratospheric  warming  frequency  due  to  use  of  stochastic  or  deterministic 
nonorographic  gravity- wave  drag  (Charlton  et  al.  2007).  Nonetheless,  Figure  8  highlights  the 
potential  for  the  stochastic  scheme  to  generate  different  GCM  behavior  than  the  determin¬ 
istic  version  by  more  random  forcing  that  can  seed  large  irreversible  changes  via  nonlinear 
interactions  and  feedbacks.  The  differences  in  the  summer  hemisphere  in  Figure  8b,  which 
are  also  larger  than  those  in  the  summer  hemisphere  in  Figure  6d,  are  probably  due  to 
interhemispheric  coupling  through  a  modified  mesospheric  pole-to-polc  residual  circulation 
caused  by  modified  gravity-wave  driving  in  northern  winter  due  to  the  differently  disturbed 
winter  stratospheres  in  each  simulation  (Becker  and  Fritts  2006). 

Variability  in  the  deterministic  scheme’s  drag  comes  solely  from  variability  in  resolved 
GCM  winds  and  temperatures,  which  in  turn  peaks  at  planetary  wavenumbers.  Thus, 
the  deterministically  parameterized  nonorographic  gravity-wave  drag  forces  variability  most 
strongly  at  the  gravest  GCM  wavenumbers.  By  contrast,  the  stochastic  scheme’s  drag  also 
varies  significantly  and  randomly  from  point  to  point  in  both  space  and  time  (Figure  4), 
and  thus  could  force  significant  variability  at  the  smallest  space-time  scales  of  the  GCM. 
Since  these  smallest  GCM  scales  can  be  unreliable  and  can  alias  to  larger  scales  (Lander 
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and  Hoskins  1997),  such  small-scale  forcing  might  have  undesirable  side  effects  on  the  GCM 
simulations. 

To  investigate  this,  in  each  nature  run  the  GCM’s  instantaneous  global  spectral  fields 
were  saved  at  0000  UTC  on  every  model  day.  Figure  9a  plots  mean  total  kinetic  energy 
spectra  at  0.055  hPa  for  June  (2007-2009)  for  the  stochastic  and  deterministic  runs,  along 
with  the  divergent  and  vortical  contributions,  as  computed  directly  from  these  daily  spherical 
harmonic  spectral  coefficients  of  GCM  vorticity  and  divergence  (e.g.,  Koshyk  et  ah  1999). 
In  computing  these  long-term  means,  we  also  computed  standard  deviations  at  each  total 
wavenumber,  which  are  plotted  in  Figure  9b.  In  the  row  beneath,  we  plot  the  ratio  of  these 
spectral  distributions  between  the  GCM  fields  using  stochastic  and  deterministic  drag,  for 
both  the  mean  spectra  (Figure  9c)  and  their  standard  deviations  (Figure  9d).  Overall,  we 
do  not  see  any  noticeable  change  in  either  the  shape  or  intensity  of  the  mean  GCM  kinetic 
energy  spectra  between  the  stochastic  and  deterministic  simulations.  The  only  clear  signal 
is  in  Figure  9d,  which  shows  that  the  standard  deviation  of  the  divergent  component  of 
the  GCM  kinetic  energy  is  systematically  larger  at  the  highest  total  wavenumbers  in  the 
simulation  with  stochastic  gravity-wave  drag. 

To  study  these  trends  as  a  function  of  altitude,  we  averaged  the  spectral  means  and 
standard  deviations  at  each  height  over  the  band  of  total  wavenumbers  60-75  (shaded  in 
Figure  9d).  Figure  10  plots  the  height  variation  of  the  ratio  of  the  mean  spectra  and  their 
standard  deviations  between  the  stochastic  and  deterministic  nature  runs.  Once  again,  the 
only  clear  trend  with  altitude  is  a  slight  overall  increase  in  the  standard  deviation  of  divergent 
kinetic  energy  in  the  GCM  at  high  total  wavenumbers  in  the  upper  stratosphere  and  meso¬ 
sphere  when  using  the  stochastic  scheme,  although  mean  divergent  kinetic  energy  averaged 


16 


over  height  is  also  slightly  higher  systematically  in  Figure  10a.  Why  a  similar  enhancement 
in  the  standard  deviation  of  small-scale  GCM  vorticity  is  not  evident  in  Figure  10b  is  not 
clear.  It  may  be  that  this  stochastic  gravity-wave  forcing  produces  regional  imbalances  that 
lead  to  the  rapid  regeneration  of  secondary  resolved  gravity  waves  (divergent  energy)  in  the 
GCM  that  propagate  away  to  restore  a  balanced  vorticity  field  (e.g.,  Vadas  et  al.  2003).  In 
any  event,  these  results  indicate  that  highly  stochastic  and  intermittent  parameterizations 
of  gravity-wave  forcing,  in  addition  to  producing  reliable  GCM  climate,  do  not  appear  to  be 
accompanied  by  any  major  unrealistic  increases  in  small-scale  dynamical  variability  within 
the  GCM. 

6.  Discussion 

While  there  have  been  occasional  efforts  to  parameterize  nonorographic  gravity-wave  drag 
stochastically  in  GCMs  (e.g.,  Dunkerton  1982a;  Piani  et  al.  2004),  the  parameterizations 
currently  used  in  production  GCM  configurations  are  exclusively  multi-wave  deterministic 
formulations.  A  strong  practical  motivation  for  the  current  stochastic  approach  to  param¬ 
eterizing  nonorographic  gravity-wave  drag  is  to  reduce  the  computational  expense  of  these 
multi-wave  deterministic  schemes  in  GCMs.  When  the  65-wave  WACCM  3.0  scheme  of 
Garcia  et  al.  (2007)  was  implemented  in  NOGAPS-ALPHA,  for  instance,  it  alone  consumed 
between  20-35%  of  the  forecast  model’s  total  run  time.  The  extension  of  NOGAPS-ALPHA 
into  the  middle  atmosphere  involves  additional  model  layers  and  physics  that  are  slated  for 
future  transition  to  the  operational  NOGAPS.  As  a  numerical  weather  prediction  (NWP) 
prototype,  new  NOGAPS-ALPHA  features  compete  for  scarce  computational  resources  with 
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many  other  potential  upgrades  to  NOGAPS  with  potentially  greater  immediate  impact  on 
NWP  (e.g.,  higher  horizontal  resolution,  more  tropospheric  observations  for  assimilation, 
etc.).  Thus,  to  be  a  viable  candidate  for  near-term  transition  to  operations,  parameteriza- 
tions  of  nonorographic  gravity-wave  drag  must  be  both  accurate  and  computationally  cheap 
relative  to  the  total  run  time  of  the  system. 

The  single-wave  stochastic  analogue  of  the  65-wave  scheme  developed  here  should  ide¬ 
ally  yield  close  to  a  65-fold  increase  in  computational  speed.  As  discussed  in  Appendix  A, 
it  was  implemented  here  within  the  existing  parameterization  code,  which  contains  signif¬ 
icant  additional  overhead  associated  with  internal  calculations  of  different  meteorological 
profiles  and  time-mean  statistics,  and  thus  a  65-fold  speed  increase  cannot  be  expected. 
Nonetheless,  without  any  additional  effort  to  further  optimize  this  parameterization  subrou¬ 
tine,  the  single-wave  stochastic  option  yields  an  order  of  magnitude  increase  in  the  speed 
of  this  subroutine  relative  to  the  deterministic  65-wave  version  in  the  NOGAPS-ALPHA 
nature  runs  reported  here.  The  single-wave  stochastic  scheme  now  consumes  between  1-4% 
of  the  forecast  model’s  total  run  time  and  thus  becomes  a  viable  transition  candidate.  The 
computational  expense  issue  is  unique  neither  to  this  particular  parameterization  nor  to  this 
particular  GCM.  For  example,  in  their  implementation  of  the  deterministic  nonorographic 
gravity-wave  drag  scheme  of  Scinocca  (2003)  in  the  European  Centre  for  Medium-Range 
Weather  Forecasts  Integrated  Forecast  System  (ECMWF  IFS),  Orr  et  al.  (2010)  discretized 
their  source  momentum-flux  spectrum  using  ngw  =  80  individual  waves.  Due  to  the  resultant 
computational  expense,  they  found  it  necessary  to  update  the  tendency  from  this  scheme 
every  two  hours  only,  to  reduce  the  overall  computational  burden  to  ~3%  of  the  total  run 
time.  Our  stochastic  approach  allows  us  to  achieve  similar  or  greater  computational  savings 
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without  resorting  to  a  reduced  space-time  physics  grid  in  the  GCM. 

A  convenient  aspect  of  the  stochastic  formulation  developed  here  is  the  close  relation¬ 
ship  that  is  retained  to  the  antecedent  deterministic  schemes,  given  that  the  latter  schemes 
are  now  common  in  GCMs  and  have  been  exhaustively  tuned  over  many  years  to  yield 
realistic  middle  atmosphere  climate.  Using  simple  conversion  relations,  our  offline  single- 
column  simulations  showed  that  identical  time-mean  mean-flow  accelerations,  heating  rates, 
and  diffusivities  could  be  generated  using  a  straightforward  stochastic  analogue  of  a  tuned 
determi  nisi  tic  scheme.  More  importantly,  when  implemented  within  a  GCM,  nearly  indistin¬ 
guishable  zonal-mean  drag  and  zonal-mean  climate  were  produced  in  July,  for  example.  Such 
reproducible  GCM  climate  responses  were  not  assured  given  the  potentially  large  nonlinear 
feedbacks  involved  in  transitioning  from  drag  that  is  smooth  and  deterministic,  to  drag  that 
is  noisy  and  random  on  small  space-time  scales.  These  findings  essentially  accord  with  those 
of  McLandress  and  Scinocca  (2005),  who  found  that  GCMs  were  remarkably  insensitive  to 
the  precise  ways  in  which  nonorographic  gravity-wave  momentum  fluxes  were  deposited  as 
a  function  of  height  in  different  deterministic  schemes. 

The  stochastic  scheme  produces  random  drag  variability  at  the  smallest  space-time  scales 
of  the  GCM  that  is  entirely  absent  in  the  deterministic  parent  scheme.  In  essence,  this 
variability  now  makes  explicit  in  the  GCM  the  inherent  gravity-wave  intermittency  that 
is  parameterized  in  the  deterministic  scheme  using  the  bulk  scaling  factor  e.  A  body  of 
literature  has  highlighted  potential  advantages  of  an  explicit  representation  of  such  random 
intrinsic  intermittency  in  GCM  parameterizations  (Palmer  2001;  Palmer  et  al.  2005;  Wilks 
2008),  such  as  more  realistic  ensemble  spread  and  variability  (Buizza  et  al.  1999,  2005; 
Teixeira  and  Reynolds  2008;  Reynolds  et  al.  2008)  and  mean  error  reduction  via  more  realistic 
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population  of  different  climate  and  weather  regimes  (Molteni  and  Tibaldi  1990;  Jung  et  al. 
2005).  Figure  8,  for  example,  showed  very  different  mean  winter  polar  stratospheric  winds  in 
January  between  the  stochastic  and  deterministic  nature  runs  due  to  stratospheric  warmings 
in  the  stochastic  GCM  simulation  that  did  not  occur  in  the  deterministic  simulation.  These 
2.5  year  runs  are  far  too  short  to  deduce  any  systematic  differences  in  stratospheric  warming 
frequency.  Nonetheless,  these  results  are  consistent  with  simple  conceptual  models  that  show 
how  small  amounts  of  random  gravity  wave  forcing  can  trigger  large  regime  transitions  that 
generate  stratospheric  warmings  that  do  not  occur  in  a  corresponding  deterministic  model 
without  random  gravity- wave  forcing  (Birner  and  Williams  2008). 

Despite  intense  forcing  at  the  smallest  space-time  scales  of  the  GCM  by  the  stochastic 
gravity-wave  drag  scheme,  our  GCM  simulations  did  not  reveal  any  large  systematic  in¬ 
creases  in  mean  kinetic  energy  at  the  smallest  GCM  scales  relative  to  the  corresponding 
deterministic  simulation,  although  a  small  systematic  increase  in  the  standard  deviation  of 
divergent  kinetic  energy  at  small  resolved  scales  was  clearly  evident  (Figure  9).  This  finding 
may  explain  why  GCMs  that  numerically  suppress  realistic  kinetic  energy  at  small  space- 
time  scales  see  greatest  improvements  not  via  stochastic  parameterization  alone,  but  also 
by  explicitly  injecting  additional  stochastic  kinetic  energy  back  into  these  smallest  resolved 
GCM  scales  (Jung  et  al.  2005;  Berner  et  al.  2009;  Charron  et  al.  2010). 

Recent  work  has  also  found  that  mesoscale  GCM  kinetic  energy  in  the  mesosphere  is 
intrinsically  chaotic  and  stochastic  due  to  the  dominance  of  divergent  gravity  wave  motions 
that  have  fast  decorrelation  times  (Liu  et  al.  2009;  Nezlin  et  al.  2009).  Thus,  the  random 
stochastic  forcing  of  the  smallest  mesospheric  space-time  scales  of  a  GCM  using  a  stochastic 
nonorographic  gravity-wave  drag  parameterization  may  in  fact  mimic  the  true  stochastic 
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nature  of  these  dynamics  in  the  GCM. 


7.  Conclusions 

We  have  presented  a  simple  methodology  for  generating  an  explicitly  stochastic  analogue 
of  an  existing  deterministic  multi-wave  parameterization  of  nonorographic  gravity-wave  drag 
that  should  be  easily  applicable  to  other  GCM  parameterizations  of  gravity-wave  drag.  Our 
approach  maintains  a  close  association  to  the  original  deterministic  scheme,  such  that  the 
stochastic  version  is  implemented  here  as  an  option  and  minor  modification  of  the  original 
deterministic  parameterization  code.  Through  the  use  of  simple  scaling  terms,  we  show  how 
the  stochastic  analogue  reproduces  identical  time-mean  drag,  diffusion  and  heating  rates  to 
the  deterministic  parent  scheme,  which  greatly  simplifies  replacing  the  latter  with  the  former 
in  GCMs  using  the  existing  tuned  parameter  settings  of  the  deterministic  antecedent. 

When  implemented  in  a  GCM,  our  single-wave  stochastic  analogue  of  the  65-wave  WACCM 
3.0  deterministic  nonorographic  gravity-wave  drag  scheme  produced  largely  identical  zonal- 
mean  climate  and  very  similar  spectral  energy  distributions  in  long-term  nature  runs  to 
those  from  corresponding  runs  using  the  original  deterministic  scheme.  In  addition  to  repro¬ 
ducing  very  similar  GCM  climate  and  variability,  the  stochastic  scheme  yields  the  following 
additional  beneficial  features: 

•  an  order-of-magntitude  reduction  in  computational  expense; 

•  explicit  parameterization  of  gravity-wave  intermittency,  which  largely  replaces  the 
tuned  bulk  intermittency  factor  e  in  the  deterministic  scheme; 
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•  stochastic  GCM  variability  that  can  realistically  increase  ensemble  spread  and  reduce 
mean  climate  biases. 

Given  these  benefits  along  with  no  apparent  disadvantages  to  date,  we  are  now  routinely 
parameterizing  nonorographic  gravity-wave  drag  in  NOGAPS-ALPHA  using  the  single-wave 
stochastic  parameterization  outlined  here.  Having  made  this  change,  new  parameterization 
possibilities  now  open  up.  For  example,  there  are  emerging  observations  of  gravity-wave 
momentum-flux  intermittency  (e.g.,  Hertzog  et  al.  2008),  which  could  now  be  used  to  con¬ 
strain  the  explicit,  stochastic  variability  of  the  parameterization  more  realistically  (see,  e.g., 
Figure  5).  Similarly,  there  are  other  parameters  in  the  scheme  besides  wave  phase  speeds 
that  could  be  converted  from  deterministic  to  stochastic  variables.  Obvious  candidates  are 
those  that  are  poorly  constrained  observationally  or  theoretically  and  which  are  likely  to 
vary  quasi-randomly  rather  than  having  set  values,  such  as  the  background  momentum  flux 
rfc*  and  the  launch  pressure  level  psrc,  among  others. 

A  longer-term  goal  is  to  transition  from  the  crude  background  sources  used  here  to  more 
physical  parameterizations  of  nonorographic  gravity-wave  drag  from  specific  tropospheric 
sources  inferred  from  GCM  fields,  such  as  deep  convection  and  jet  instabilities  (Charron 
and  Manzini  2002;  Beres  et  al.  2005;  Richter  et  al.  2010).  One  might  assume  that  physical 
source  models  would  naturally  lead  back  to  a  deterministic  parameterization  approach,  but 
that  may  be  unlikely.  GCM  parameterizations  of  convection  remain  challenging  and  are 
themselves  moving  towards  explicitly  stochastic  models  (Plant  and  Craig  2008;  Teixeira 
and  Reynolds  2008),  which  would  naturally  require  a  stochastic  parameterization  of  the 
drag  due  to  parameterized  gravity  waves  emanating  from  such  stochastic  convective  sources 
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in  a  GCM.  Similarly,  the  jet  instabilities  and  circulations  that  generate  gravity  waves  are 
highly  resolution  dependent  (e.g.,  Scinocca  and  Ford  2000;  Plougonven  and  Snyder  2007). 
When  these  instabilities  are  properly  resolved  by  the  GCM,  they  likely  generate  resolved 
waves  that  do  not  need  to  be  parameterized  (e.g.,  O’Sullivan  and  Dunkerton  1995).  Thus, 
the  inadequately  resolved  nonorographic  sources  of  unresolved  gravity-wave  momentum  flux 
that  require  parameterization  in  GCMs  may  continue  to  be  more  usefully  parameterized 
stochastically  to  reflect  this  inherent  uncertainty  and  variability.  Finally,  it  should  be  noted 
that  aspects  of  orographic  gravity-wave  drag  also  appear  to  be  inherently  stochastic  (Doyle 
and  Reynolds  2008;  Eckermann  et  al.  2010),  which  may  motivate  explicitly  stochastic  schemes 
to  replace  the  deterministic  parameterizations  of  orographic  gravity-wave  and  flow-blocking 
drag  that  are  currently  used  in  GCMs  (e.g.,  Palmer  2001). 
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APPENDIX  A 


Nonorographic  Gravity  Wave  Drag  Module 

In  2006,  NASA,  NCAR  and  NRL  began  a  collaborative  project  to  develop  gravity-wave- 
drag  parameterization  jointly  for  all  three  member  GCMs:  the  Goddard  Earth  Observing 
System  Version  5  (GEOS5),  WACCM  and  NOGAPS-ALPHA,  respectively.  The  scheme 
described  by  Kiehl  et  al.  (1996)  formed  the  basis  for  parameterized  nonorographic  gravity 
wave  drag  initially  implemented  within  both  WACCM  and  GEOS5.  That  common  parame¬ 
terization  subsequently  diverged.  NCAR,  inter  alia,  added  more  parameterized  waves,  used 
a  wider  phase-speed  distribution  and  changed  the  source-level  momentum  flux  for  WACCM 
(Garcia  et  al.  2007).  NASA  adopted  different  source  functions  and  modifed  the  propagation 
and  dissipation  modules  for  use  in  GEOS5.  The  Garcia  et  al.  (2007)  formulation  in  WACCM 
was  later  implemented  in  NOGAPS-ALPHA,  where  it  too  was  modified  and  tuned  for  data 
assimilation  applications  (see  section  3  of  Eckermann  et  al.  2009).  At  the  same  time,  the 
version  in  WACCM  underwent  further  large  independent  changes  (Richter  et  al.  2010). 

It  soon  proved  impossible  for  each  center  to  continually  integrate  into  its  version  of  the 
code  all  the  new  features  emerging  at  the  other  two  centers,  especially  as  the  codes  at  each 
center  became  more  dissimilar  with  time.  This  spurred  a  programming  effort  to  combine  the 
three  different  versions  of  the  code  at  each  center  into  a  single  common  parameterization  that 
all  insititutions  could  then  implement  in  their  GCMs  and  develop  jointly  from  a  common  code 
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and  repository.  This  so-called  “team  scheme”  was  carefully  coded  from  scratch  to  adhere 
rigidly  to  the  11  “plug  compatibility”  rules  proposed  by  Kalnay  et  al.  (1989)  to  facilitate 
easier  exchange  of  parameterizations  among  modeling  centers. 

The  biggest  change  relative  to  the  antecedent  codes  was  the  creation  of  separate  set-up 
and  running  subroutines  (Rule  2),  the  former  entirely  new.  This  new  set-up  subroutine  is 
called  just  once  at  the  start  of  a  GCM  run.  Through  two  simple  input  labels  -  a  “model” 
and  an  “experiment”  identifier  -  a  series  of  specific  statements  and  parameter  settings  are 
activated  that  define  subsequent  behavior  of  the  gravity-wave-drag  subroutine.  The  “model” 
label  identifies  a  particular  GCM  by  activating  gravity-wave  drag  options  and  parameter 
values  used  in  that  GCM,  and  deactivating  all  other  features  used  in  other  GCMs.  The 
“experiment”  label  activates  a  secondary  series  of  settings  that  activate  preprogrammed 
“tuned”  parameter  values  for  a  particular  GCM  experiment,  configuration  or  resolution. 

Given  that  the  “team  scheme”  aims  to  be  a  community  multi-institution  GCM  resource, 
the  new  set-up  subroutine  offers  many  advantages.  For  example: 

•  backwards  compatibility:  older  code  and/or  tuned  parameter  settings  can  be  retained 
and  easily  reactivated  to  rerun  historical  GCM  experiments  or  configurations. 

•  faster  transitions:  new  physics  options  developed  for  one  center’s  GCM  now  become 
immediately  available  for  other  centers  to  activate  and  test  in  their  GCM. 

•  greater  flexibility:  all  physics  options  are  available  and  can  now  be  easily  mixed  and 
matched,  or  simply  deactivated. 

•  separation  of  “tuning”  from  code  development:  casual  GCM  users  seeking  only  to 
“tune”  the  parameterization  in  a  GCM  need  only  to  make  a  few  simple  edits  to  the 
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set-up  subroutine.  The  core  physics  subroutines  of  the  parameterization  should  now 
never  need  to  be  edited  by  anyone  other  than  parameterization  developers. 

While  this  integrated  capability  comes  with  an  inevitable  increase  in  the  overall  length 
and  complexity  of  the  code  relative  to  it’s  antecedents,  this  is  more  than  compensated  for 
by  these  and  other  advantages. 

Offline  and  online  tests  of  the  team  scheme  have  verified  exact  reproduction  of  the  results 
of  the  three  antecedent  nonorographic  gravity- wave-drag  codes  previously  run  at  NCAR, 
NASA  and  NRL,  so  that  each  center  can  now  use  the  team  scheme  without  any  change  in 
the  tuned  gravity-wave  drag  settings  they  have  always  used  in  their  GCM.  There  is  also  no 
significant  speed  penalty  of  the  new  code  relative  to  those  antecdent  codes  in  timing  tests 
to  date  using  the  NOGAPS- ALPHA  GCM. 

This  team  scheme  also  integrates  the  different  orographic  gravity  wave  drag  schemes 
used  at  each  center,  specifically:  (a)  an  orographic  gravity-wave  scheme  based  on  McFarlane 
(1987)  used  by  NCAR  and  NASA  (Kiehl  et  al.  1996),  and  (b)  a  Webster  et  al.  (2003)  param¬ 
eterization  of  orographic  gravity-wave  and  flow-blocking  drag  used  in  NOGAPS-ALPHA.  In 
the  team-scheme,  source-level  orographic  gravity-wave  fluxes  from  either  scheme  are  sent  to 
the  same  common  propagation  and  dissipation  modules  used  in  the  nonorographic  gravity- 
wave  calculations. 

The  stochastic  parameterization  of  nonorographic  gravity  wave  drag  outlined  in  this 
paper  has  also  been  implemented  in  the  team  scheme  as  a  new  option,  and  was  used  to 
generate  all  the  offline  and  online  (GCM)  results  presented  in  this  paper. 
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Fig.  1.  Background  momentum  flux  (mPa)  based  on  two  different  functions  F((j),t )  in 
(2)  used  by  (a)  WACCM  3.0  (rfe*  =  7  mPa)  and  (b)  NOGAPS-ALPHA  (rfe*  =  10  mPa). 
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(b)  Stochastic  Momentum  Flux  Distribution 

/  \  ♦V=1  ; 

\  c  =  RAN  DOM  [-80  m  s  \+80  m  s 

^  \c  =  RANDOM[-40  m  s  \+40  m  s'1]  “ 


Fig.  2.  Tsrc(cj)/Tb  for  (a)  ngw  =  65  discretization  of  Garcia  et  al.  (2007)  and  ngw  =  9 
discretization  of  Kiehl  et  al.  (1996),  and  (b)  the  corresponding  stochastic  analogues. 
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Fig.  3.  Vertical  profiles  of  (a)  zonal  and  (b)  meridional  winds  input  to  gravity  wave  drag 
scheme,  which  returns  (c)  zonal  and  (d)  meridional  mean-flow  accelerations,  (e)  vertical 
diffusion  coefficients,  and  (f)  heating  rates.  As  labeled  in  panel  (c),  the  gray  curves  show 
results  from  the  WACCM  3.0  scheme,  while  the  remaining  curves  show  output  from  the 
stochastic  analogue  averaged  over  1,  10,  100,  1000,  and  10000  separate  calls. 
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Fig.  4.  (a)  zonal  and  (b)  meridional  mean-flow  accelerations  from  the  WACOM  3.0  gravity 
wave  drag  scheme  (gray  curve)  and  the  10000-point  mean  from  the  stochastic  analogue  (black 
curve)  with  corresponding  standard  deviations  shown  as  error  bars. 
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Fig.  5.  As  in  Figure  4a,  but  now  showing  means  and  standard  deviations  from  the  stochastic 
analogue  for  nsgw  values  of  (a)  1  (same  as  Figure  4a),  (b)  9,  (c)  21,  (d)  81,  (e)  201,  and  (f) 
801.  In  each  case  the  number  of  calls  was  adjusted  to  yield  a  total  of  10000  waves. 
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(a)  Control  (no  MGWD)  July  (2007-2009)  (b)  Deterministic  MGWD  July  (2007-2009) 


Fig.  6.  Zonal-mean  zonal  winds  (m  s_1)  for  July  2007-2009  from  NOGAPS-ALPHA  nature 
runs:  (a)  control  run  without  nonorographic  gravity  wave  drag,  and  runs  using  the  (b) 
deterministic  and  (c)  stochastic  parameterizations  of  nonorographic  gravity  wave  drag  with 
equivalent  settings  to  produce  the  same  time-mean  drag.  Differences  in  the  zonal-mean  zonal 
winds  between  the  stochastic  and  deterministic  simulations  are  plotted  in  (d). 
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(a)  Deterministic  3U/3t:  July  (2007-2009) 


(b)  Stochastic  3U/3t:  July  (2007-2009) 
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(b)  Stochastic  3T/3t:  July  (2007-2009) 
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Fig.  7.  Zonal-mean  zonal  mean-flow  accelerations  (top  row,  nr  s-1  day-1)  and  heating  rates 
(bottom  row,  K  day-1)  due  to  parameterized  nonorographic  gravity  wave  drag  averaged  for 
July  2007-2009  from  NOGAPS-ALPHA  nature  runs  using  (left  column)  deterministic  and 
(middle  column)  stochastic  parameterizations.  Differences  between  the  two  parameteriza- 
tions  are  plotted  in  the  right  column  panels. 
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Fig.  8.  Zonal-mean  zonal  winds  (m  s_1)  for  January  2008-2009  from  NOGAPS-ALPHA 
nature  runs:  (a)  control  run  without  nonorographic  gravity  wave  drag,  and  runs  using  the 
(b)  deterministic  and  (c)  stochastic  parameterizations  of  nonorographic  gravity  wave  drag 
with  equivalent  settings  to  produce  the  same  time-mean  drag.  Differences  in  the  zonal-mean 
zonal  winds  between  the  stochastic  and  deterministic  simulations  are  plotted  in  (d). 
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(c)  Mean  Energy  Spectral  Ratios 


(d)  Standard  Deviation  Ratios 


Fig.  9.  (a)  mean  kinetic  energy  spectra  and  (b)  its  standard  deviation  at  each  total 

wavenumber,  at  0.055  hPa,  averaged  from  daily  0000  UTC  spectral  NOGAPS-ALPHA  GCM 
fields  for  June  2007-2009  for  nature  runs  in  which  nonorographic  gravity- wave  drag  was  pa¬ 
rameterized  deterministically  (gray  curves)  and  stochastically  (black  curves).  The  vortical 
(dashed  curves)  and  divergent  (dotted  curves)  contributions  to  total  kinetic  energy  (solid 
curves)  are  also  shown.  Panels  below  show  the  ratios  of  the  stochastic  to  the  deterministic 
spectral  curves  in  the  panels  above.  Panel  (d)  shades  the  total  wavenumber  range  60-75 
used  to  form  profile  means  in  Figure  10. 
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(a)  Mean  Spectral  Ratios 
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(b)  Standard  Deviation  Ratios 
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Fig.  10.  Ratios  of  total  (solid),  divergent  (dotted)  and  vortical  (dashed)  contributions  to 
(a)  mean  kinetic  energy  and  (b)  its  standard  deviation  in  the  NOGAPS-ALPHA  GCM  fields 
for  June  2007-2009  using  parameterized  stochastic  and  deterministic  nonorographic  gravity- 
wave  drag.  Each  energy  value  is  averaged  over  the  total  wavenumber  range  60-75,  with  the 
resulting  ratios  between  the  stochastic  and  deterministic  GCM  fields  plotted  as  a  function 
of  height.  Values  greater  (less)  than  unity  imply  increased  (decreased)  kinetic  energy  at 
wavenumbers  60-75  in  the  GCM  fields  with  stochastic  nonorographic  gravity-wave  drag, 
relative  to  those  with  deterministic  nonorographic  gravity-wave  drag. 
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